TD Dynamique adaptative - MDet

Modèle écologique

On considère un modèle de compétition Dans notre modèle, 21 phénotypes dont le trait X est distribué sur un intervalle allant de 0 à 2, sont en compétition. Chaque population de chaque phénotype possède un taux de croissance intrinsèque identique r. La densité maximale est spécifique à chaque phénotype et dépend uniquement de la valeur de son trait selon la formule :

\[ K(x) = max(K_0 - λ(x-x_0)^2, 0) \] L’intensité de la compétition entre deux phénotypes i et j vaut : \[ a(x_i, x_j) =exp(-\frac{(x_i-xj)^2}{2σ^2}) \]

La dynamique de chaque population s’écrit alors :

\[ \frac{dn_i}{dt} = r.n_i.(\sum_{j=1}^{n} \frac{a(x_i, x_j).n_j}{K(x_i)} ) \] Pour l’ensemble des manipulations, on fixe certaines valeurs par défaut aux paramètres.

#Parametres
i = 21 # Nombre d'especes
Xmin = 0
Xmax = 2
X0 = (Xmax - Xmin)/2
K0 = 1
lambda = 1
sigma = 1 #etalement de la competition
r = 1

On implémente les fonctions mathématiques du modèle.

############ Fonctions #############

K_cap = function(x, K0_ = K0, lambda_= lambda, x0 = X0){
   (K0_ - lambda_*(x-x0)**2)
}

a = function(x1,x2, sigma_a = sigma){ 
  exp( (-1/2*(x2-x1)**2)/sigma_a**2 )
}

fitness = function(x1, x2, r_ = r ){
  r_*(1-a(x1, x2)*(K_cap(x1)/K_cap(x2)))
}

Fitness d’invasion

Afin de comprendre ce qu’implique notre modèle avec les valeurs de ses paramètres initiaux, nous représentons graphiquement la valeur de la fitness relative dans le cas de l’apparition d’un nouveau phénotype. Ces figures en trois dimensions fournissent le paysage adaptatif de notre modèle.

## Set range and domain of plot
x1_interval  <- seq(0.1, 1.9, length.out = 25);
x2_interval  <- seq(0.1, 1.9, length.out = 25);

## Interpolate surface

z  <- outer(x1_interval,x2_interval,
            FUN = fitness)

p  <- persp(x1_interval,x2_interval,z, theta = 30, phi = 20,
            col = "lightblue", shade = 0.4, ticktype = "detailed")

plot_ly() %>% add_surface(x = ~x1_interval, y = ~x2_interval, z = ~z)

Le Pairwise Invasibily Plot (PIP) fournit un résumé de la situation en représentant seulement le signe de la différence des fitness.

## Set range and domain of plot
x1_interval  <- seq(from = 0.1, to = 1.9, length.out = 500);
x2_interval  <- seq(from = 0.1, to =1.9, length.out = 500);

## Interpolate surface

z  <- outer(x1_interval,x2_interval,
            FUN = fitness)

couleurs <- ifelse(z > 0, "lightblue", "black")

matrice_couleurs = couleurs
df <- data.frame(x = rep(x1_interval, each = ncol(matrice_couleurs)),
                 y = rep(x2_interval, times = nrow(matrice_couleurs)),
                 couleur = as.vector(matrice_couleurs))

ggplot(df, aes(x = x, y = y, fill = couleur)) +
  geom_tile() +
  scale_fill_identity() +
  labs(title = "PIP", xlab = 'x1', ylab = 'x2') +
  theme_minimal()+
  labs(fill = "Légende\nNoir = Négatif\nBleu = Positif")

On déduit de ce PIP que le trait x = 1 est instable et attractif, il remplit les conditions pour être un ‘branching point’ ce qui sera visualisé dans la suite (partie sur la dynamique évolutive avec simulation de mutations).

Simulations Lotka-Volterra

Pour simuler les compétitions, le modèle est implémenté sous forme matricielle. La dynamique est simulée pour un temps de 1000 et pour deux valeurs de sigma.

############### Simuler LV pour un nombre arbitraire d'espece ############

time_max = 1000

#Conditions initiales
x = seq(from = 0, to = 2, length.out = i) # Valeurs des phenotypes
k = c()
for (i in 1:length(x)){
  k = c(c(k),c(K_cap(x[i])))
}# K de chaque espece

n = rep(1/i,i) #densite de pop initiale

# On met sous forme de matrice
N0 = t(as.matrix(n))
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)
X = matrix(rep(x,i), ncol = i, nrow = i, byrow = TRUE)

#Calcul de la matrice M
D = t(X) - X


#Resolution du systeme d'ODE
model <- function(t,N,sigma){
  
  M = exp(-0.5*D^2/(sigma^2))/t(K)
  
  dN <- r * N * (1 - N %*% M)
  
  return(list(dN))
}

#ode <-ode(N0,0:time_max,model,parms = sigma)

Evolution de la densité des populations selon la valeur du trait pour sigma = 5 (faible étalement de la compétition)

sigma = 5
# Representation

ode <-ode(N0,0:time_max,model,parms = sigma)

ode <- as.data.frame(ode)
#colnames(ode) <- c("t",'x1','x2','x3')
ode_plot = ode%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ggplot(ode_plot)+
  geom_line(aes(time, dens, col = sp_ID))+
  labs(col = "Population \nidentity")+
  xlab("Time") + ylab("Population density")+
  ggtitle("σ = 5")

Etalement de la compétition avec sigma = 5

ode_plot$sp_ID <- as.numeric(ode_plot$sp_ID)
ode_plot$trait <- x[as.numeric(ode_plot$sp_ID)]

ggplot(ode_plot)+
  geom_raster(aes(trait, time,  fill=dens**2))+
  scale_fill_gradient2(low = "white" ,
                       high = "red")+
  labs(fill = "Population \ndensity")+
  xlab("Trait") + ylab("Time")+
  ggtitle("σ = 5")

Evolution de la densité des populations selon la valeur du trait pour sigma = 0.3 (fort étalement de la compétition)

sigma = 0.3
# Representation

ode <-ode(N0,0:time_max,model,parms = sigma)

ode <- as.data.frame(ode)
#colnames(ode) <- c("t",'x1','x2','x3')
ode_plot = ode%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ggplot(ode_plot)+
  geom_line(aes(time, dens, col = sp_ID))+
  labs(col = "Population \nidentity")+
  xlab("Time") + ylab("Population density")+
  ggtitle("σ = 0.3")

Etalement de la compétition avec sigma = 0.3

ode_plot$sp_ID <- as.numeric(ode_plot$sp_ID)
ode_plot$trait <- x[as.numeric(ode_plot$sp_ID)]

ggplot(ode_plot)+
  geom_raster(aes(trait, time,  fill=dens**2))+
  scale_fill_gradient2(low = "white" ,
                       high = "red")+
  labs(fill = "Population \ndensity")+
  xlab("Trait") + ylab("Time")+
  ggtitle("σ = 0.3")

On observe que le nombre de phénotypes qui persistent dans le temps dépend notamment de la valeur de sigma. Pour des faibles valeurs, un plus grand nombre de phénotypes se maintiennent.

Le script et le graphique suivant permettent d’étudier le nombre de phénotypes se maintenant en fonction de la valeur de sigma. Les valeurs de sigma testées sont comprises entre 0.1 et 2.

seuil <- 0.001
sigma_bank = seq(0.1,2,by = 0.05)

nb_especes <- c()

for (sigma_i in sigma_bank){{
  
  ode2 <-ode(N0,0:time_max,model,parms = sigma_i)
  nb = sum(ode2[nrow((ode2)),-1] > seuil)
  }
  nb_especes <- c(c(nb_especes),c(nb))
  
}

res = tibble(sigma_bank,nb_especes)

ggplot()+
  geom_point(data = res, aes(x = sigma_bank, y = nb_especes))+
  labs(title = " ",
       x = "Sigma",
       y = "Nombre d'espèces au régime stationnaire")

Modélisation de mutations : dynamique adaptative

Dans cette partie, on modélise la dynamique adaptative sur le temps long (T) d’une population en considérant que p = 10% de la population mute à intervalle de temps régulier (t). La population initiale est monomorphe et l’étalement de la compétition (sigma = 0.2) est important, on s’attend donc à ce que sur le temps long beaucoup de phénotypes apparaissent et se maintiennent

#Implementation des mutations

#parametres
p = 0.1 #taux de la pop qui mute
Temps = 10000 #temps de la simulation
t = 400 #pas de temps entre chaque mutation
n = rep(0,i) 
n[4] <- 1 # densites de pop initiale (une seule espece)
N0 <- t(as.matrix(n))
sigma = 0.2

#Initialisation des valeures de K
k = c()
for (i in 1:length(x)){
  k = c(c(k),c(K_cap(x[i])))
}
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)


A = exp(-0.5*D^2/(sigma^2))
M = A/t(K)



# Evolution sans mutation sur le premier pas de temps t
ode <-ode(N0,1:t,model,parms = sigma, method = 'lsoda')
ode_mutation <- ode


seuil = 0.05 #densité min on considère qu'un trait avec une densité < seuil a une densité de 0


for (z in 2:(Temps/t)){
  
  #Elimination des traits sous le seuil
  tail <- as.matrix(tail(ode,1)[,1:i+1])
  for (j in 1:length(tail)){
    if(tail[j] < seuil){
      tail[j] <- 0
    }
  }
  N0 <- t(tail)
  
  # Mutation des traits restants
  for (index in which(tail != 0)){
    a = runif(1)
    if (a<0.5){ #mutation vers la gauche
      N0[max(0,index-1)] <- p*N0[index]
      N0[index] <- (1-p)*N0[index]
    } else { #mutation vers la droite
      N0[min(i,index+1)] <- p*N0[index]
      N0[index] <- (1-p)*N0[index]
    }
   
  }
  
  
  ode <-ode(N0,(z*t):(z*t+t),model,parms = sigma, method='lsoda')
  
  ode_mutation <- rbind(ode_mutation,ode) 

}

ode_mutation_plot = as.data.frame(ode_mutation)%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ode_mutation_plot$sp_ID <- as.numeric(ode_mutation_plot$sp_ID)
ode_mutation_plot$trait <- x[as.numeric(ode_mutation_plot$sp_ID)]

ggplot(ode_mutation_plot)+
   geom_tile(aes(trait, time, fill = dens))+
   scale_fill_gradient2(low = "white" ,
                        high = "red") +
  labs(fill = "Population \ndensity")+
  xlab("Trait") + ylab("Time")+
  ggtitle("Adaptative dynamics")

On constate effectivement en premier lieu une convergence vers la valeur de trait optimale = 1 (Attractivité de x = 1), puis des branchements évolutif vers d’autres valeures de phénotypes (Instabilité de x = 1).